function positive_definite_basics()
% POSITIVE_DEFINITE_BASICS 正定和半正定矩阵基础
% 
% 介绍正定和半正定矩阵的基本概念、性质和判定方法

fprintf('=== 正定和半正定矩阵基础 ===\n\n');

%% 1. 基本定义
fprintf('1. 基本定义\n');
fprintf('对于实对称矩阵A:\n');
fprintf('- 正定: x^T*A*x > 0 对所有非零向量x\n');
fprintf('- 半正定: x^T*A*x ≥ 0 对所有向量x\n');
fprintf('- 负定: x^T*A*x < 0 对所有非零向量x\n');
fprintf('- 半负定: x^T*A*x ≤ 0 对所有向量x\n\n');

% 示例矩阵
A_pos = [2, 1; 1, 2];  % 正定
A_semi = [1, 1; 1, 1];  % 半正定
A_neg = [-2, -1; -1, -2];  % 负定
A_indef = [1, 2; 2, 1];  % 不定

matrices = {A_pos, A_semi, A_neg, A_indef};
names = {'正定', '半正定', '负定', '不定'};

for i = 1:length(matrices)
    A = matrices{i};
    fprintf('%s矩阵:\n', names{i});
    disp(A);
    
    % 计算特征值
    eigenvals = eig(A);
    fprintf('特征值: [%.4f, %.4f]\n', eigenvals);
    
    % 判定正定性
    if all(eigenvals > 0)
        type = '正定';
    elseif all(eigenvals >= 0)
        type = '半正定';
    elseif all(eigenvals < 0)
        type = '负定';
    elseif all(eigenvals <= 0)
        type = '半负定';
    else
        type = '不定';
    end
    fprintf('判定结果: %s\n\n', type);
end

%% 2. 正定性的等价条件
fprintf('2. 正定性的等价条件\n');
fprintf('对于实对称矩阵A，以下条件等价:\n');
fprintf('(1) A是正定的\n');
fprintf('(2) 所有特征值为正\n');
fprintf('(3) 所有主子式为正\n');
fprintf('(4) 存在可逆矩阵P使得A = P^T*P\n');
fprintf('(5) Cholesky分解存在\n\n');

% 测试矩阵
A_test = [4, 2, 1; 2, 3, 0.5; 1, 0.5, 2];
fprintf('测试矩阵 A =\n');
disp(A_test);

% 条件1：定义检查（通过特征值）
eigenvals_test = eig(A_test);
fprintf('条件1 - 特征值: [%.4f, %.4f, %.4f]\n', eigenvals_test);
is_positive_def = all(eigenvals_test > 0);
fprintf('所有特征值为正: %s\n', mat2str(is_positive_def));

% 条件2：主子式检查
principal_minors = zeros(3, 1);
for k = 1:3
    principal_minors(k) = det(A_test(1:k, 1:k));
end
fprintf('主子式: [%.4f, %.4f, %.4f]\n', principal_minors);
all_positive_minors = all(principal_minors > 0);
fprintf('所有主子式为正: %s\n', mat2str(all_positive_minors));

% 条件3：Cholesky分解
try
    L_chol = chol(A_test, 'lower');
    fprintf('Cholesky分解成功\n');
    fprintf('L =\n');
    disp(L_chol);
    
    % 验证 A = L*L^T
    reconstruction_error = norm(L_chol * L_chol' - A_test, 'fro');
    fprintf('重构误差: %e\n', reconstruction_error);
    
catch ME
    fprintf('Cholesky分解失败: %s\n', ME.message);
end

%% 3. 半正定矩阵的性质
fprintf('\n3. 半正定矩阵的性质\n');

% 半正定矩阵示例
B = [1, 1, 1; 1, 1, 1; 1, 1, 1];
fprintf('半正定矩阵 B =\n');
disp(B);

eigenvals_B = eig(B);
fprintf('特征值: [%.4f, %.4f, %.4f]\n', eigenvals_B);
fprintf('秩: %d\n', rank(B));

% 半正定矩阵的平方根
[V_B, D_B] = eig(B);
D_sqrt = diag(sqrt(max(diag(D_B), 0)));  % 处理数值误差
B_sqrt = V_B * D_sqrt * V_B';

fprintf('矩阵平方根 B^(1/2) =\n');
disp(B_sqrt);

% 验证 (B^(1/2))^2 = B
sqrt_error = norm(B_sqrt * B_sqrt - B, 'fro');
fprintf('平方根验证误差: %e\n', sqrt_error);

%% 4. 正定矩阵的运算性质
fprintf('\n4. 正定矩阵的运算性质\n');

A1 = [3, 1; 1, 2];
A2 = [2, 0.5; 0.5, 3];

fprintf('正定矩阵 A1 =\n'); disp(A1);
fprintf('正定矩阵 A2 =\n'); disp(A2);

% 性质1：正定矩阵的和仍是正定的
A_sum = A1 + A2;
fprintf('A1 + A2 =\n'); disp(A_sum);
eigenvals_sum = eig(A_sum);
fprintf('和的特征值: [%.4f, %.4f]\n', eigenvals_sum);
fprintf('和是正定的: %s\n', mat2str(all(eigenvals_sum > 0)));

% 性质2：正定矩阵的逆仍是正定的
A1_inv = inv(A1);
fprintf('A1的逆 =\n'); disp(A1_inv);
eigenvals_inv = eig(A1_inv);
fprintf('逆的特征值: [%.4f, %.4f]\n', eigenvals_inv);
fprintf('逆是正定的: %s\n', mat2str(all(eigenvals_inv > 0)));

% 性质3：合同变换保持正定性
P = [1, 2; 0, 1];  % 可逆矩阵
A_congruent = P' * A1 * P;
fprintf('合同变换 P^T*A1*P =\n'); disp(A_congruent);
eigenvals_cong = eig(A_congruent);
fprintf('合同变换后的特征值: [%.4f, %.4f]\n', eigenvals_cong);
fprintf('合同变换后是正定的: %s\n', mat2str(all(eigenvals_cong > 0)));

%% 5. 条件数和数值稳定性
fprintf('\n5. 条件数和数值稳定性\n');

% 良态正定矩阵
A_good = diag([1, 2, 3, 4]);
cond_good = cond(A_good);
fprintf('良态正定矩阵条件数: %.2f\n', cond_good);

% 病态正定矩阵
epsilon = 1e-6;
A_ill = diag([1, 2, 3, epsilon]);
cond_ill = cond(A_ill);
fprintf('病态正定矩阵条件数: %.2e\n', cond_ill);

% Cholesky分解的稳定性
L_good = chol(A_good, 'lower');
L_ill = chol(A_ill, 'lower');

fprintf('良态矩阵Cholesky因子条件数: %.2f\n', cond(L_good));
fprintf('病态矩阵Cholesky因子条件数: %.2e\n', cond(L_ill));

%% 6. 正定矩阵的构造
fprintf('\n6. 正定矩阵的构造\n');

% 方法1：A^T*A形式
n = 4;
X = randn(n, n);
A_construct1 = X' * X;
fprintf('方法1 - A^T*A构造:\n');
fprintf('最小特征值: %.6f\n', min(eig(A_construct1)));

% 方法2：对角占优
A_construct2 = randn(n, n);
A_construct2 = A_construct2 + A_construct2';  % 对称化
A_construct2 = A_construct2 + (abs(min(eig(A_construct2))) + 1) * eye(n);  % 确保正定
fprintf('方法2 - 对角占优构造:\n');
fprintf('最小特征值: %.6f\n', min(eig(A_construct2)));

% 方法3：指定特征值
desired_eigenvals = [1, 2, 3, 4];
V_random = orth(randn(n, n));  % 随机正交矩阵
A_construct3 = V_random * diag(desired_eigenvals) * V_random';
fprintf('方法3 - 指定特征值构造:\n');
actual_eigenvals = sort(eig(A_construct3));
fprintf('期望特征值: [%.1f, %.1f, %.1f, %.1f]\n', desired_eigenvals);
fprintf('实际特征值: [%.6f, %.6f, %.6f, %.6f]\n', actual_eigenvals);

%% 7. 应用：二次型优化
fprintf('\n7. 应用：二次型优化\n');
fprintf('最小化二次函数 f(x) = (1/2)*x^T*A*x - b^T*x\n\n');

% 优化问题参数
A_opt = [2, 1; 1, 3];
b_opt = [1; 2];

fprintf('Hessian矩阵 A =\n'); disp(A_opt);
fprintf('线性项 b = [%.1f; %.1f]\n', b_opt);

% 检查凸性
eigenvals_opt = eig(A_opt);
fprintf('Hessian特征值: [%.4f, %.4f]\n', eigenvals_opt);
is_convex = all(eigenvals_opt > 0);
fprintf('函数是凸的: %s\n', mat2str(is_convex));

if is_convex
    % 解析解
    x_optimal = A_opt \ b_opt;
    f_optimal = -0.5 * b_opt' * x_optimal;
    
    fprintf('最优解: x* = [%.6f; %.6f]\n', x_optimal);
    fprintf('最优值: f* = %.6f\n', f_optimal);
    
    % 验证最优性条件：梯度为零
    gradient = A_opt * x_optimal - b_opt;
    fprintf('最优点梯度: [%.2e; %.2e]\n', gradient);
end

%% 8. 矩阵不等式
fprintf('\n8. 矩阵不等式\n');
fprintf('Löwner偏序: A ⪰ B 当且仅当 A - B 半正定\n\n');

A_ineq1 = [3, 1; 1, 2];
A_ineq2 = [2, 0.5; 0.5, 1];

fprintf('矩阵 A =\n'); disp(A_ineq1);
fprintf('矩阵 B =\n'); disp(A_ineq2);

diff_matrix = A_ineq1 - A_ineq2;
fprintf('A - B =\n'); disp(diff_matrix);

eigenvals_diff = eig(diff_matrix);
fprintf('A - B的特征值: [%.4f, %.4f]\n', eigenvals_diff);

if all(eigenvals_diff >= -1e-10)  % 考虑数值误差
    fprintf('A ⪰ B (A大于等于B)\n');
else
    fprintf('A ⪯ B (A小于等于B)\n');
end

%% 9. 正定矩阵的几何解释
fprintf('\n9. 正定矩阵的几何解释\n');
fprintf('正定矩阵定义椭球面\n\n');

% 2D椭球示例
A_ellipse = [2, 0.5; 0.5, 1];
fprintf('椭球矩阵 A =\n'); disp(A_ellipse);

% 特征值分解得到椭球的主轴
[V_ellipse, D_ellipse] = eig(A_ellipse);
eigenvals_ellipse = diag(D_ellipse);
semi_axes = 1 ./ sqrt(eigenvals_ellipse);

fprintf('主轴方向 (特征向量):\n');
disp(V_ellipse);
fprintf('半轴长度: [%.4f, %.4f]\n', semi_axes);

% 绘制椭球
theta = linspace(0, 2*pi, 100);
unit_circle = [cos(theta); sin(theta)];

% 变换到椭球
A_sqrt_inv = V_ellipse * diag(1./sqrt(eigenvals_ellipse)) * V_ellipse';
ellipse_points = A_sqrt_inv * unit_circle;

figure('Name', '正定矩阵的几何解释');
plot(ellipse_points(1, :), ellipse_points(2, :), 'b-', 'LineWidth', 2);
hold on;
plot(unit_circle(1, :), unit_circle(2, :), 'r--', 'LineWidth', 1);

% 绘制主轴
scale = max(semi_axes);
for i = 1:2
    axis_vec = scale * semi_axes(i) * V_ellipse(:, i);
    plot([-axis_vec(1), axis_vec(1)], [-axis_vec(2), axis_vec(2)], ...
         'g-', 'LineWidth', 2);
end

axis equal;
grid on;
xlabel('x_1');
ylabel('x_2');
title('椭球 x^T A x = 1');
legend('椭球', '单位圆', '主轴', 'Location', 'best');

%% 10. 数值检验工具
fprintf('\n10. 数值检验工具\n');

% 创建检验函数
function result = check_positive_definite(A, method)
    switch method
        case 'eigenvalues'
            eigenvals = eig(A);
            result = all(eigenvals > 1e-12);
        case 'cholesky'
            try
                chol(A);
                result = true;
            catch
                result = false;
            end
        case 'minors'
            n = size(A, 1);
            result = true;
            for k = 1:n
                if det(A(1:k, 1:k)) <= 1e-12
                    result = false;
                    break;
                end
            end
        case 'quadratic'
            % 随机测试二次型
            n_tests = 100;
            n = size(A, 1);
            result = true;
            for i = 1:n_tests
                x = randn(n, 1);
                if x' * A * x <= 1e-12
                    result = false;
                    break;
                end
            end
    end
end

% 测试不同的检验方法
test_matrices = {[2, 1; 1, 2], [1, 1; 1, 1], [2, 3; 3, 2]};
test_names = {'正定', '半正定', '不定'};
methods = {'eigenvalues', 'cholesky', 'minors', 'quadratic'};

fprintf('正定性检验方法比较:\n');
fprintf('矩阵类型    特征值    Cholesky    主子式    二次型\n');
fprintf('--------    ------    --------    ------    ------\n');

for i = 1:length(test_matrices)
    A_check = test_matrices{i};
    fprintf('%-8s    ', test_names{i});
    
    for j = 1:length(methods)
        result = check_positive_definite(A_check, methods{j});
        fprintf('%-8s  ', mat2str(result));
    end
    fprintf('\n');
end

end